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Abstract 


The relationship between variations in the Earth’s orbit and the periodic expansion and re¬ 
treat of continental ice sheets over the last 2.5 million years (2.5myr) remains poorly de¬ 
fined. Spectral analysis of globally-averaged benthic <5 18 0 isotope ratios, a proxy record of 
global ice volume, yields peaks at frequencies that align with variations in each of eccen¬ 
tricity, precession and obliquity. This is interpreted as strong evidence that orbital variability 
is the key driver of large-scale climate change over the Pleistocene epoch. However, geo¬ 
logical evidence shows ice age cycles to be globally synchronised, with both hemispheres 
experiencing glacial or interglacial conditions simultaneously. This is inconsistent with the 
structure of orbital insolation forcing, which we show to be highly asymmetric between the 
hemispheres. This asymmetry is produced principally through variations in the precession 
parameter, which dominates insolation variability in almost all latitudes and seasons. Obliq¬ 
uity is relevant in the very high latitudes, and interacts with precession to produce strong 
regional insolation anomalies in these regions. Global annual mean insolation variability is a 
function of the eccentricity parameter alone, but of negligible influence to the climate system 
as a whole. The more significant role played by eccentricity is in modulating the amplitude 
of precessional forcing. Using the Globally Resolved Energy Balance (GREB) model, this 
thesis presents a regional, seasonal and global annual mean deconstruction of the Earth’s 
climate response to orbital variability. We find that the oceans, atmospheric circulation, the 
ice-albedo feedback and the water vapour feedback act to strongly amplify and seasonally 
rearrange orbital insolation forcing. Atmospheric CO 2 variability is shown to influence the 
time-scale behaviour of the global annual mean climate evolution, but is not the determina¬ 
tive forcing mechanism. Using these insights, we introduce a simple insolation model that 
better defines the relationship between the orbital parameters and the global annual mean 
response of the GREB model. This model, based on the obliquity-and-eccentricity modula- 
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tion of precession (OEMP), scales and combines variations in each of the three parameters. 
It demonstrates that the Southern Hemisphere dominates the global mean response of the 
GREB climate. This is the inverse of the observed glacial signal, which is shown to be 
driven by variations in insolation over the Northern Hemisphere. For all configurations of 
the GREB model presented in this thesis, globally synchronised ice ages are not successfully 
simulated, suggesting that additional climate feedbacks are required to reproduce this key 
feature of glacial-interglacial cycles. Ice sheet dynamics and an interactive global carbon 
cycle, both of which are not simulated in the GREB model, are likely to be highly relevant. 

Keywords: Ice age cycles, glacial-interglacial, orbital forcing, Milankovitch cycles, insola¬ 
tion, climate model. 
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Introduction 


The Earth’s paleoclimate history over the Pleistocene epoch (last 2.5myr) has been dom¬ 
inated by semi-regular cycles of large-scale glacial expansion and retreat. During glacial 
conditions, continental ice sheets extended over much of northern Eurasia and North Amer¬ 
ica, and Antarctic temperatures were up to 9°C colder than today (Petit et al., 1999). Benthic 
5 18 0 ratios from foraminifera in marine sediment cores are sensitive to global ice sheet ex¬ 
tent and provide clear physical evidence of this ice age signal (Bintanja, van de Wal, & 
Oerlemans, 2005). Spectral analysis of the globally averaged <5 18 0 record yields peaks con¬ 
sistent with the key periodicities of three parameters that describe the Earth’s orbital geome¬ 
try (Hays, Imbrie, & Shackleton, 1976; Shackleton, 2000). These variations in eccentricity, 
precession and obliquity strongly influence the regional and seasonal distribution of inso¬ 
lation that reaches the top of the atmosphere (Milankovitch, 1941). Despite the clear link 
between orbital variability and large-scale glaciation over the Pleistocene, the exact nature 
of this relationship remains poorly understood (Huybers, 2011; Willed, Ganopolski, Calov, 
Robinson, & Maslin, 2015; Elkibbi & Rial, 2001). A key area of uncertainty is the in¬ 
consistency between the observed glacial signal, which is globally synchronised, and the 
dominant structure of the orbital forcing, which produces a forcing anomaly of opposite sign 
in each hemisphere (Mercer, 1984; Denton et al., 1999). The prevailing theory holds that 
summer (usually measured at the June solstice) insolation at 65°N ultimately governs large- 
scale glacial variability, as this describes the latitude and season at which the sensitivity of 
Northern Hemisphere ice sheets is maximised (Oerlemans, 1991; Deblonde, Peltier, & Hyde, 
1992). However, the relationship between 65°N summer insolation and the <5 18 0 signal is 
poor, as demonstrated in figure 0.1. Although interglacial periods (peaks in figure 0.1) often 
coincide with large insolation anomalies (e.g. blue shaded areas), large interglacials also 
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occur during periods of only moderate insolation peaks (e.g. red shaded areas). Similarly, 
insolation variability at 65°N is dominated by ~21,000-year (21kyr) fluctuations in preces¬ 
sion, whereas observed ice ages over the last lmyr vary at a period of roughly lOOkyr. This 
implies that the Earth’s climate must respond to the underlying 65°N insolation variabil¬ 
ity in a strongly non-linear fashion (Doughty et al., 2015). Although climate models have 
been effective in reproducing individual features of Pleistocene glacial variability (e.g. Pol¬ 
lard, 1983; Paillard, 1998; Tziperman, Raymo, Huybers, & Wunsch, 2006; Huybers, 2006; 
Parrenin & Paillard, 2012; Tzedakis et al., 2012), none have been successful in establish¬ 
ing a comprehensive relationship between the orbital forcing signal and the Earth’s climate 
response (Willeit, Ganopolski, Calov, & Brovkin, 2019). 
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Figure 0.1: a Globally-averaged benthic t) 18 0 stack for the past 800kyr. b Daily mean insolation at 65°N on 
June 21. Marine Isotope sub-Stages displayed above a. Arrows denote interglacials, defined as the point at 
which the benthic 3 ls O crosses a 3.68%o threshold. Modified from Tzedakis et al. (2017). 


Eccentricity, precession and obliquity, collectively known as the Milankovitch mecha¬ 
nisms, describe the features of the Earth’s orbit that influence the seasonal and spatial distri¬ 
bution of insolation reaching top of the Earth’s atmosphere. Figure 0.2 illustrates the physical 
property recorded by each parameter. Eccentricity defines the shape of the Earth’s orbit, with 
larger values indicating a more elliptical path around the sun. This means eccentricity also 
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defines the difference between perihelion (the point in the Earth’s orbit at which it is closest 
to the sun) and aphelion (the point at which the Earth-Sun distance is largest). This is rel¬ 
evant because the entire Earth will receive relatively more solar radiation at the perihelion 
of its orbit than at aphelion. The precession parameter records the ‘wobbling’ motion of 
the Earth’s axis, and therefore defines how the seasons are positioned relative to perihelion 
and aphelion. Under the current arrangement, the Northern Hemisphere winter solstice falls 
near to perihelion (see figure 0.2). The Earth travels faster at this point of its orbit, meaning 
that the Northern Hemisphere winter is both warmer and shorter than the long-term average. 
Precession is measured in degrees, with a value of 180°equivalent to the statement ‘perihe¬ 
lion occurs exactly on the March equinox’. A value of 270° therefore refers to perihelion 
coinciding with the June solstice; the September equinox for 360° and the December solstice 
for a value of 90°. The obliquity parameter measures the tilt of the Earth’s axis relative to its 
orbital path, which varies between 22.08°and 24.48°. A more pronounced tilt will enhance 
the seasonal cycle for all latitudes, with the largest effects over the very high latitudes. 



Figure 0.2: A visual representation of the the eccentricity (a), precession (b) and obliquity (c) orbital param¬ 
eters. Eccentricity describes the shape of the Earth’s orbit around the sun, which varies from more elliptical 
(solid line) for large values, to more circular (dashed line) for small values. Precession records how the seasons 
are positioned relative to perihelion and aphelion. Currently, the Northern Hemisphere winter solstice (Decem¬ 
ber 21/22) falls near perihelion. The time taken (in days) for the Earth to travel between each orbital position is 
also indicated. Obliquity provides a measure of the tilt of the Earth’s axis of rotation, which has a current value 
of 23°27'. Modified from Berger et al. (2007). 


This thesis analyses the relationship between orbital insolation variability and the ob¬ 
served ice age cycles of the Pleistocene epoch. In particular, it discusses the inconsistency 
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between the dominant asymmetric insolation forcing pattern and the globally synchronised 
climate response preserved in the paleoclimate record. Chapter 1 focuses on the regional and 
seasonal structure of Milankovitch forcing, and will deconstruct the overall insolation pattern 
into the relative contribution of each orbital parameter. A series of sensitivity experiments 
using the GREB model will then be presented in chapter 2, addressing the key processes 
relevant to the overall response of the Earth system. The insights gained from these first two 
chapters are brought together in chapter 3, which introduces a simple insolation model to 
discuss how the three orbital parameters interact to produce a global mean response in both 
the GREB climate and in the <5 18 0 record. 


Methods 

The numerical solutions to variations in eccentricity, precession and obliquity over the past 
5myr are taken from Huybers and Eisenman (2006). The combined effect of this orbital 
variability in moderating the latitudinal and seasonal distribution of insolation reaching the 
top-of-the-atmosphere is calculated according to the formula in Berger (1978, appendix). 
Individual parameter patterns are obtained by calculating the insolation anomaly produced 
when each parameter is artificially held at a constant value over a specified period. The 
Globally-Resolved Energy Balance (GREB) model is a simple climate model consisting of 
three layers: land and ocean surface, atmosphere and sub-surface ocean. The GREB model 
divides the Earth surface into a 3.75°x 3.75“latitude-longitude grid and simulates physical 
processes with a numerical timestep of 12 hours. The main processes simulated by the 
GREB model are incoming solar (shortwave) and outgoing thermal (longwave) radiation, at¬ 
mospheric heat transport by isotropic diffusion and advection, a hydrological cycle and heat 
uptake by the subsurface ocean. The sub-surface ocean extends to three times the maximum 
depth of the seasonally-shifting mixed-layer (climatology from Lorbacher, Dommenget, 
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Niiler, & Kohl, 2006). The GREB model boundary conditions include a land-sea mask, 
topography and cloud cover climatology from Rossow and Schiffer (1991). 

Prior to the analysis of this thesis, the GREB model had not been applied to paleocli- 
mate scenarios. The original model used to perform simulations was a prototype code, and 
significant testing and development efforts were required to produce accurate simulations on 
paleoclimate timescales. A successful implementation of the GREB paleoclimate scenarios 
incorporates prescribed variations in the three orbital parameters taken from Huybers and 
Eisenman (2006). The simplicity of the GREB model is such that global climate simula¬ 
tions run at a speed of 100,000 years per day on a laptop computer. This is considerably 
faster than the majority of paleoclimate models. The GREB model also allows for the de- 
construction of the overall climate system by turning OFF individual processes for specific 
sensitivity experiments. Two key feedback mechanisms simulated in the GREB model are 
relevant on glacial (>10kyr) timescales: the water vapour feedback and the ice albedo feed¬ 
back. This latter mechanism is parameterised from the surface temperature over the land 
and ocean (Dommenget & Floter, 2011). For land areas, points below -10°C area considered 
completely ice-covered and points above 0°C completely ice-free, with a linear relationship 
in the transition zone. For ocean areas, these threshold temperatures are -7°C and -1.7°C. 
Sea ice variability includes associated changes in surface heat capacity, which is calculated 
as equivalent to that of the open ocean for ice free areas and to a 2m water column for ice 
covered regions. The five GREB sensitivity experiments analysed in this thesis are presented 
in table 0.1. 

The first four sensitivity experiments incorporate a constant atmospheric CO 2 concen¬ 
tration of 280ppm (pre-industrial) as a boundary condition. In addition, the All ON GREB 
experiment is repeated with prescribed atmospheric CO 2 variability from Antarctic ice core 
data (see figure 0.2). As this ice core record extends only 800kyr into the past, the data for the 
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Table 0.1: The GREB model sensitivity experiments. 


Experiment Name & 

Length 

External Forcing 

Description 

No. 




exp-1. Local 

lmyr 

Orbital 

Thermal radiation response only. Consists only of a topography-free land 

Radiation Balance 


(prescribed) 

mask and a uniformly 50m-deep ocean. No other processes considered. 

exp-2. Circulation 

lmyr 

Orbital 

All climate processes other than the two feedbacks considered. Heat 



(prescribed) 

transport in the atmosphere and heat exchange with the subsurface ocean. 

exp-3. Ice-Albedo 

lmyr 

Orbital 

As for exp-2, but with the ice albedo feedback. The heat capacity 



(prescribed) 

influence of sea ice is simulated. 

exp-4. Water Vapour 

lmyr 

Orbital 

As for exp-3, but with the water vapour feedback. All processes and 

(All ON) 


(prescribed) 

feedbacks considered (All ON). 

exp-5. GREB-CO 2 

lmyr 

Orbital & CO 2 

(prescribed) 

As for exp-4, but with prescribed atmospheric CO 2 concentrations based 

on ice core data from Liithi et al. (2008). 


first 200kyr of the GREB-CO 2 run is a duplicate of the data from 731.7kyr BP to 531.7kyr 
BP. This range was selected to minimise the discontinuity at the transition between the ob¬ 
served and reproduced data. As the first 200kyr of the GREB-CCb simulation is included 
only for the model spin-up and is not directly analysed, the non-physical CO 2 evolution over 
this period does not affect the accuracy of the model results. 



Age (kyrs) 


Figure 0.2: Prescribed atmospheric CO 2 variability in the GREB-CO 2 sensitivity experiment. Note: CO 2 
concentration from lmyr BP - 800kyr BP is a duplicate of years 731.7kyr BP - 531.7kyr BP. Data from (Liithi 
et al., 2008) 




Chapter 1 


Orbital Insolation Forcing 


The close link between 5 18 0 isotope ratios in marine sediment cores and the three orbital 
parameters strongly suggests a direct causal relationship between the two. However, the 
global annual mean influence of orbital variability on the amount of solar radiation reaching 
the top of the atmosphere is minimal. This chapter therefore provides a comprehensive 
analysis of the regional, seasonal and time-scale behaviour of orbital insolation variability. 
This will provide the basis for better understanding how the climate interacts with this forcing 
to produce a global annual mean signal, as preserved in the geological record. 

1.1 Orbital Timescales 

Eccentricity, precession and obliquity each describe particular features of the Earth’s orbital 
configuration, and vary on millennial time scales due to the gravitational influence of nearby 
planets. Numerical solutions of the configuration of these parameters provide a reliable 
record for the past several million years (Berger, 1978). Figure 1.1 presents these numerical 
solutions over the past lmyr, as well as the key periods of variability for each parameter. The 
eccentricity parameter varies between a minimum of 0.0046 and a maximum of 0.057, with 
dominant periods of 413kyr and 102kyr. Physically, precession describes the calendar day 
(and therefore season) at which perihelion occurs. As this is measured as an angle, plotting 
the time evolution of this parameter is of more physical relevance if the sine transformation 
of the raw data is taken. This removes discontinuities that occur as the threshold of 360°is 
crossed. The physical interpretation of this arbitrary transformation is that a value of 1 
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Figure 1.1: Variations in the three orbital parameters over the past lmyr. Precession shown as the sine transfor¬ 
mation of the raw data. Obliquity measured in degrees; eccentricity and precession values are dimensionless. 
The log-linear power spectrum for each parameter is also shown, with dominant periods of variability labelled. 
Data from Huybers and Eisenman (2006). 

and -1 describes perihelion coinciding with the December and June solstices, respectively. 
Precession is characterised by spectral power at two main frequencies: 23kyr and 19kyr. The 
third parameter, obliquity, is dominated by cycles with a ~41kyr periodicity. 

1.2 Orbital Insolation Variability 

Variations in the three orbital parameters influence the global, regional and seasonal distri¬ 
bution of insolation that reaches the top of the Earth’s atmosphere. Figure 1.2 presents the 
overall strength of this orbital insolation variability. It shows that variability for all extra- 
tropical latitudes is essentially a function of summer insolation, which is roughly three times 
greater than mean winter variability. This means that orbital variability in general acts to ei¬ 
ther enhance or dampen the seasonal cycle. The effects are therefore maximised in the high 
latitudes where the seasonal cycle is strongest. Variability is maximised in the very high 
latitudes in summer, with a peak standard deviation of ~25 W/m 2 . 
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Figure 1.2: Daily mean insolation variability due to variations in eccentricity, precession and obliquity over the 
past lmyr. Data from Huybers and Eisenman (2006). 


1.3 Patterns of Variability: The Three Parameters 

The overall insolation response to orbital forcing can be deconstructed into the patterns of 
variability produced by each of the parameters. Figure 1.3 presents the mean insolation 
anomaly pattern produced by a maximum, minimum and mean value for eccentricity and 
obliquity. Precession is defined as an angle and has no maximum, minimum or mean value. 
These key values are therefore arbitrarily defined as a June perihelion, a June aphelion and 
a March equinox, respectively. It is important to note that this analysis shows the mean 
signal over a million-year period, and that the insolation anomaly for any one year will not 
correspond exactly with any of the patterns presented in figure 1.3. 

1.3.1 Eccentricity 

Global annual mean insolation variability is a function of the eccentricity parameter alone 
(see figure 1.3). Although the precise eccentricity anomaly pattern for any one year will de¬ 
pend on the position of perihelion, in mean terms a large eccentricity will produce a positive 
insolation anomaly (figure 1.3a) and a small eccentricity will produce a negative anomaly 
(figure 1.3d). However, the magnitude of this forcing anomaly is negligible: less than 2 
W/m 2 for all possible configurations. The seasonal and latitudinal effects are similarly weak, 
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Figure 1.3: The mean insolation anomalies produced by each of the three orbital parameters (columns 1-3) 
at three key configurations over the past lmyr. Row 1: maximum value; row 2: minimum value; row 3: 
mean value. Precession configurations are defined as a June perihelion, a June aphelion and a March equinox, 
respectively. Global annual mean anomaly also shown. Note: colourbar for eccentricity differs from precession 
and obliquity. Data from Huybers and Eisenman (2006). 


with a peak anomaly in the summer high latitudes of ~3 W/m 2 . The direct effect of eccen¬ 
tricity alone therefore plays a minor role in contributing to overall insolation variability over 
the Pleistocene. 


1.3.2 Precession 

Figure 1.3 presents an important property of the precession parameter: a June perihelion con¬ 
figuration is identical to a December aphelion configuration. Thus, the insolation anomaly 
pattern for a June perihelion (figure 1.3b) is the inverse of the June aphelion (figure 1.3e) 
anomaly. In general, there are two key configurations of the precession parameter: perihelion 
falling on an equinox (hereafter equinox precession forcing) or a solstice (hereafter solstice 
precession forcing). All precessional values describe one of these two arrangements or the 
transition between them. Unlike eccentricity, the influence of precession is strictly seasonal 
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in nature, with no net global annual mean insolation anomaly for all possible values (see 
figure 1.3). However, the mean latitudinal and seasonal effects are the strongest of all three 
parameters, peaking at 30 W/m 2 . As perihelion must occur each year, the mean influence 
of precession is equally strong for all possible values that it may take. This characteristic is 
not shared by eccentricity or obliquity, for which insolation variability is maximised during 
extremes, but weak for intermediate values. Precession is also the only parameter to have 
any significant effect on insolation outside of the high latitudes. A key feature of the pre¬ 
cession signal is that it is in-phase between hemispheres but out-of-phase between seasons. 
This means that for any given calendar day, the precession anomaly will be of the same sign 
in both hemispheres. However, any one latitude will receive a forcing anomaly of opposite 
sign during opposite seasons (or exactly six months later). 

The pattern of equinox precession forcing (e.g. figure 1.3h) shows a peak anomaly on 
the equator and is equally strong and of opposite sign on either hemisphere. On an equinox, 
the strongest effects will occur directly over the equator and decrease with increasing lati¬ 
tude. This structure is repeated six months later at the opposing equinox, with the sign of 
the forcing reversed. The forcing pattern will therefore be exactly mirrored between the 
hemispheres for each individual calendar day, meaning that averaged over an entire year, 
the insolation forcing for all latitudes is zero. This is only the case when perihelion (and 
therefore aphelion) occurs exactly on one of the two equinoxes. 

Under solstice precession configuration (e.g. figure 1.3b), perihelion occurs when the 
sun is positioned directly overhead either the Tropic of Cancer or the Tropic of Capricorn. 
Thus, unlike the equinox precessional configuration, a solstice precession forcing results in 
a distinct asymmetry in the strength of the insolation forcing between the hemispheres for 
each individual calendar day. This means that the insolation anomaly over both hemispheres 
will peak in summer, with a positive (negative) anomaly for the hemisphere that experiences 
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perihelion (aphelion) in this season. This also means that the forcing anomaly experienced by 
the winter hemisphere will be much smaller in magnitude, enhancing the seasonal insolation 
cycle for all latitudes. For all extra-tropical latitudes, this has the effect of creating a net 
non-zero annual mean insolation anomaly. As the seasonal insolation cycle becomes more 
significant with increasing latitude, the magnitude of this annual mean insolation forcing 
also increases towards the poles, with the strongest effects occurring at latitudes higher than 
40°N/S. Thus, for all precessional configurations except for an equinox precession forcing, 
most latitudes will receive a non-zero annual mean insolation anomaly. When perihelion 
coincides exactly with one of the two solstices, this annual mean anomaly is maximised. 
However, as the forcing will be mirrored exactly between the two hemispheres, the global 
annual mean forcing will remain zero for all precessional configurations. 

1.3.3 Obliquity 

Figure 1.3c & f present the mean insolation anomaly pattern produced by a maximum and 
minimum obliquity, respectively. They show that a more pronounced axial tilt will strongly 
enhance the insolation reaching the high latitudes in summer. The opposite will occur in 
the winter hemisphere, but as the polar latitudes receive no insolation in this season, this 
winter forcing is of limited significance. Although the maximum influence of obliquity on 
insolation is comparable to the strength of precession forcing, this is restricted to the very 
high latitudes in summer and only for maximum or minimum values. A mean obliquity value 
(figure 1.3i) produces no significant mean insolation anomaly for any latitude or calendar 
day. The obliquity signal is out-of-phase between the hemispheres but in-phase seasonally 
- the opposite configuration to precession. This means that for any given value of obliquity, 
the hemispheres will receive an insolation anomaly of opposite sign on the same calendar 
day. However, each pair of latitudes north and south of the equator will receive the exact 
same forcing signal exactly six months apart. 
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1.4 Parameter Interactions 


This section describes how particular combinations of parameter configurations interact with 
each other to enhance or dampen the seasonal and latitudinal distribution of insolation vari¬ 
ability. As the precessional parameter produces the largest mean insolation anomaly (see 
figure 1.3) this section will focus on how eccentricity and obliquity interact with this domi¬ 
nant precession signal. 

Eccentricity combines with precession forcing through two key transformations. Firstly, 
the global annual mean eccentricity anomaly is combined with the seasonal precession forc¬ 
ing signal. This has the effect of shifting the precessional insolation anomaly at each latitude 
and calendar day ‘up’ (more positive/less negative) for a large eccentricity and ‘down’ (less 
positive/more negative) for small eccentricity. Secondly, the eccentricity of the Earth’s orbit 
modulates the amplitude of the precession anomaly. This is because the difference between 
perihelion and aphelion increases for a more elliptical orbital path around the sun. The com¬ 
bined influence of these two interactions is summarised in figure 1.4 (left), which plots the 
zonal annual mean insolation anomaly when perihelion occurs on the June solstice. The 
mean anomaly (black) is compared with the combined anomaly during a maximum (red) 
and minimum (blue) eccentricity value. 

Under the precessional configuration shown in 1.4, the northern hemisphere expe¬ 
riences perihelion during its summer, and therefore the zonal mean insolation anomaly is 
positive (negative) for all northern (southern) latitudes. If this scenario coincides with a 
maximum eccentricity (red), the two interactions described above (i) enhance the magnitude 
of the anomaly in both hemispheres, and (ii) shift the anomaly for each latitude towards more 
positive values. This upward shift forces the positive anomaly experienced by the northern 
hemisphere to be slightly stronger than the negative anomaly over the southern hemisphere 
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Figure 1.4: Zonal annual mean insolation anomaly produced when perihelion falls on the June solstice. Black 
line: mean over the past lmyr; red line: anomaly produced when a maximum eccentricity (left) or obliquity 
(right) interacts with a June solstice precessional configuration; blue line: anomaly produced when a minimum 
eccentricity (left) or obliquity (right) interacts with a June solstice precessional configuration. 


(~+17 W/m 2 vs. ~-14 W/m 2 ). The inverse process is repeated when a June solstice perihe¬ 
lion occurs in concert with a minimum eccentricity value, yielding anomalies of ~+7 W/m 2 
and ~-7 W/m 2 in the northern and southern hemispheres, respectively. The strength of this 
modulating effect of eccentricity on precession is such that it is often of more physical rel¬ 
evance to combine the two into a single parameter, referred to as ‘eccentricity-modulated 
precession’ or EMP (see figure 3.1 for a graphical representation). This parameter measure¬ 
ment is obtained by multiplying precession by the matching eccentricity value. 


Figure 1.4 (right) shows the zonal annual mean insolation anomaly produced from the 
interaction between a June perihelion and a maximum (red) and minimum (blue) obliquity. 
The effect of obliquity is largest in the very high latitudes, while the low latitudes receive 
an anomaly that is both smaller in amplitude and opposite in sign (see figure 1.3). This is 
reflected in figure 1.4 (right), where the insolation anomaly over the high latitudes in both 
hemispheres is shifted towards more positive values for a maximum obliquity. The opposite 
occurs during a minimum obliquity. Regardless of the value taken by obliquity, the effect of 
the interaction between the two parameters in the low latitudes (< 60°N/S) is small. 
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1.5 Insolation Magnitude and Spatial Structure 


Figure 1.5 presents the magnitude of the insolation variability produced by each parameter, 
again demonstrating that variability at almost all latitudes and seasons is dominated by pre¬ 
cession. However, figure 1.5 shows the magnitude of precession forcing during a mean ec¬ 
centricity. As eccentricity modulates the strength of precession, this parameter will produce 
even larger variability when eccentricity is large. The strength of the obliquity parameter 
peaks over the polar regions in summer, where it is roughly half as strong as for precession. 
Eccentricity produces a signal that is roughly equal in magnitude in all latitudes and seasons, 
but small relative to both precession and obliquity. 
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Figure 1.5: Mean insolation variability for each of the three parameters over the last lmyr. 


Figure 1.6 depicts an EOF deconstruction of the overall insolation pattern into the key 
underlying patterns of variability, ranked according to the explained variance of each pattern 
(expressed as a percentage). This deconstruction emphasises the dominance of the preces- 
sional parameter to the overall insolation forcing pattern, with the first two EOF modes 
corresponding almost exactly with the two key precessional configurations. Together, the 
solstice forcing pattern (EOF-1) and the equinox forcing pattern (EOF-2) explain 98% of the 
overall forcing structure. The insolation pattern of obliquity (EOF-3) is comparatively lim¬ 
ited in its ability to explain the seasonal and latitudinal variations in insolation over the past 
million years. However, the obliquity pattern will project onto EOF-1, and will therefore 
contribute to the explained variance of this mode. It is important to note that the dominant 
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spatial signal produced by orbital variability (EOF-1) is asymmetrical between the hemi¬ 
spheres. By contrast, the paleoclimate signal preserved in the geological record describes 
large and globally synchronised cycles of glacial expansion and retreat. This inconsistency 
between the principal forcing signal and the observed climate response is a key source of 
uncertainty surrounding the Pleistocene ice ages, and will form the focus of the remaining 
two chapters of this thesis. 



Figure 1.6: The first three empirical orthogonal function (EOF) modes of global insolation variability over the 
past lmyr. Explained variance (%) given in the title of each mode. 
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Chapter 2 


The Climate Response 


Paleoclimate observations show synchronised temperature variability in both hemispheres 
throughout glacial-interglacial cycles. As the dominant orbital forcing pattern is fundamen¬ 
tally different to this climate signal, the Earth system must interact with solar variability 
through strong and non-linear feedback mechanisms (Ganopolski & Brovkin, 2017). This 
chapter presents several GREB model simulations of the past lmyr to discuss both the role 
of climate feedbacks in amplifying insolation variability, and the utility of the model in re¬ 
producing the observed climate variability of the Pleistocene. 

2.1 The Global Mean GREB Climate Response 

A defining characteristic of the Pleistocene glacial cycles is the distinctive ‘saw-tooth’ time 
evolution. Observations show global ice volume to gradually build over tens of thousands 
of years, before rapidly ablating within several thousand years, with a complete cycle oc¬ 
curring roughly every lOOkyr (Clark et al., 2009; Hays et al., 1976). A key challenge for 
paleoclimate models is reproducing this time-scale behaviour and dominant lOOkyr period¬ 
icity from orbital insolation forcing (Ridgwell, Watson, & Raymo, 1999; Rial, 1999; Imbrie 
et al., 1993). Figure 2.1 compares the global annual mean time evolution of the GREB model 
and the observed glacial signal. As the 8 18 record from Lisiecki and Raymo (2005) is tuned 
to the insolation signal, the GREB model is similarly tuned to the observed signal such that 
the correlation coefficient between the two is maximised. The All ON simulation (figure 2.1, 
upper) shares several features with the 8 18 record, but notably lacks the strong low-frequency 
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variability of the latter time series. In particular, variability at the lOOkyr frequency is signif¬ 


icantly less significant in the GREB simulation relative to the observed signal. These results 
suggest that an additional feedback is required to reproduce the observed ‘saw-tooth’ time 
evolution of the Pleistocene glacial-interglacial cycles. 


(5 18 0 vs. All ON Global Annual Mean Surface Temperature 




Figure 2.1: GREB model global annual mean surface temperature anomaly compared to the observed ice age 
signal from benthic 5 ls O isotope ratios. Signal are tuned such that the correlation coefficient (r) is maximised. 
As the CO 2 record from ice core measurements extends to 800kyr BR the first 200kyr of the GREB-CO 2 
simulation is not based on realistic values and therefore omitted. Note: as the relationship between (5 18 0 and 
temperature is inverse (see reversed axis), the correlation is negative. Power spectra of (5 18 ) (blue) and both All 
ON (black) and GREB-CO 2 (red) simulations are shown in logarithmic (left) and linear (right) representations. 
(5 18 0 data from Lisiecki and Raymo (2005). 


The GREB-CO 2 simulation (figure 2.1, middle) shows a much closer relationship with 
the observed 5 18 0 signal. Spectral analysis (figure 2.1c & d) show that variability at all 
frequencies is strongly enhanced relative to the All ON simulation. Global annual surface 
temperature anomalies are accordingly stronger in amplitude for the GREB-CO 2 simulation 
compared to the All ON simulation. This suggests that atmospheric CO 2 acts as a strong 
global mean feedback mechanism to orbital variability on all frequencies. With the inclu¬ 
sion of this strong positive feedback, the GREB model is able to successfully reproduce the 
characteristic ‘saw-tooth’ time evolution and dominant lOOkyr periodicity of the ice age sig- 
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nal. This suggests that the insolation forcing signal produced by orbital variability must be 
transformed by a strong internal feedback mechanism for the successful reproduction of pa- 
leoclimate observations. The insolation signal alone is insufficient. Our results suggest that 
the global carbon cycle is capable of fulfilling this purpose. Additional feedbacks not incor¬ 
porated into the GREB model are also likely to be highly significant, including the delayed 
rebound of continental landmass to ice sheet growth and retreat (Abe-Ouchi et al., 2013) and 
dust-albedo feedbacks (Ganopolski & Calov, 2011; Ellis & Palmer, 2016). 

A second key feature of the global ice age signal is the amplitude of surface temper¬ 
ature variability over glacial cycles. However, an accurate global annual mean temperature 
record on glacial timescales is problematic, principally because 5 18 0 ratios area sensitive to 
both ice sheet extent and to deep-ocean temperatures, among other local conditions (Bintanja 
et al., 2005; Imbrie et al., 1984). Difficulties in separating these competing sensitivities in¬ 
troduces significant uncertainty, particularly on global spatial scales. A more robust surface 
temperature record for the last ~420kyr has been established using ice core measurements 
from Vostok, Antarctica (Petit et al., 1999). Figure 2.2 compares the GREB annual mean 
surface temperature over the Antarctic region (defined here as 50°S - 90°S) with this recon¬ 
structed record. 


Amplitude of the GREB Response 



Figure 2.2: Amplitude of annual mean GREB surface temperature variability in the Antarctic region (50°S - 
90°S) for the All ON (black) and GREB-CO 2 (red) simulations compared with reconstructed variability from 
ice core data. Labels show the standard deviation for each time series. Data from Petit et al. (1999). 

The amplitude of the GREB model temperature response to orbital forcing is roughly 
equivalent to the ice core record, with a standard deviation of ~2.3 compared to ~2.9 for the 
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observed time series. Furthermore, there is little difference between the All ON and GREB- 


CO 2 simulations, suggesting that high latitude temperature variability is more sensitive to 
orbital fluctuations than to greenhouse forcing. Ice core temperatures are significantly more 
negative than the GREB anomalies, likely due to to presence of large ice sheets that are not 
simulated in the model experiments. Strictly in terms of amplitude, however, the GREB 
model is largely successful in reproducing the temperature variability contained in the ice 
core record. 

2.2 Deconstruction of the GREB Climate Response 

The sensitivity of the GREB climate to each of the three orbital parameters is principally con¬ 
trolled by four key processes and feedbacks: the heat capacity of the oceans, atmospheric cir¬ 
culation, the ice albedo feedback and the water vapour feedback. Using the GREB model, we 
can begin from a Local Radiation Balance sensitivity experiment, consisting only of a land 
surface mask and a uniform 50m-deep ocean, and progressively add each of these processes 
to the model to deconstruct the individual influence of each in enhancing or dampening the 
base insolation forcing. 

Figure 2.3 (left) shows the power spectra of global annual mean temperature vari¬ 
ability for each sensitivity experiment. It shows that each climate process increases global 
annual mean variability at all frequencies. In this way, an increasingly complex climate 
will strongly amplify the underlying insolation forcing due to orbital variability. Further¬ 
more, this variability is concentrated around the three main frequencies, which correspond 
to (from low-to-high frequency): eccentricity (~100kyr period), obliquity (~41kyr period) 
and the two main precessional periods (~23kyr and ~19kyr). Additional peaks at periods of 
~10.5kyr and ~7kyr are sub-harmonics of the precessional parameter, and appear because 
variability produced by this parameter is not perfectly described by a sine or cosine function. 
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As discussed in chapter 1, the insolation anomaly produced by both precession and obliquity 
is exactly zero in global annual mean. The presence of peaks at frequencies corresponding 
to these parameters therefore demonstrates that the climate must spatially and seasonally 
rearrange the direct insolation forcing to produce a signal relevant on the global scale. 


Global Annual Mean Surface Temperature Normalised by Insolation Forcing 



Figure 2.3: Power spectra for each of the GREB sensitivity experiments. Left: global annual mean surface 
temperature; right: global annual mean surface temperature normalised by insolation forcing. 

Figure 2.3 (right) shows global annual mean temperature variability normalised by 
global annual mean insolation forcing. This records the magnitude of the GREB climate 
response at each frequency, per unit forcing. Frequencies at which the climate strongly en¬ 
hances insolation variability will be larger in magnitude than frequencies at which internal 
climate processes are of reduced significance. Normalised spectral density in this represen¬ 
tation is largest for the mid-frequencies that correspond to obliquity and precessional forcing 
(see 2.3). Each successive process and feedback incorporated into the GREB model acts to 
enhance the surface temperature response to these orbital parameters. As the global annual 
mean forcing for these two parameters is essentially zero, it is on these frequencies that the 
effect of the climate is maximised. 

Interestingly, both power spectra in figure 2.3 contain peaks at the obliquity and pre¬ 
cessional frequencies for the Local Radiation Balance sensitivity experiment. This experi¬ 
ment contains only a land surface mask and a 50m-deep ocean, and so it could be expected 
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that this simplified climate system would be unable to transform a strictly seasonal forcing 
pattern (obliquity and precession) into a global annual mean climate signal (as recorded in 
figure 2.3). However, variability is strong at these frequencies, suggesting that the differing 
response of the land and the oceans to insolation forcing is highly relevant to the GREB 
model climate. 

Figure 2.4 deconstructs the surface temperature response of the Local Radiation Bal¬ 
ance simulation into the response over the land and ocean for each hemisphere in summer, 
winter and annual mean. Precession variability dominates both hemispheres in all seasons, 
producing the clear ~21 kyr variability visible in all time series. Land areas have a low 
heat capacity and respond instantaneously to this precession forcing, which is opposite in 
summer and winter. This strong seasonality over land is seen in that the surface temper¬ 
ature anomalies are exactly opposite in summer and winter. However, summer insolation 
variability is larger in amplitude and therefore dominates the annual mean response. The 
oceans, by contrast, have a much larger heat capacity and require a large forcing anomaly 
to produce a significant surface temperature response. As precessional forcing is dispropor¬ 
tionately strong in summer, ocean areas are much more sensitive to insolation variability 
during those months of the year. This signal then persists during the comparatively weak 
insolation anomaly experienced in winter. The surface temperature response over the oceans 
is therefore synchronised in both summer and winter, and in annual mean (figure 2.4b & 
d). Furthermore, the summer and winter anomalies combine to produce a far larger annual 
mean surface temperature response over the oceans relative to the land. As the majority of 
the world’s oceans are positioned below the equator, the cumulative effect of this response is 
strongest in the Southern Hemisphere (see figure 2.4, standard deviations). This integrating 
effect of the oceans allows a large annual mean variability to emerge from a strictly seasonal 
forcing signal. 
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Figure 2.4: Surface temperature variability in the GREB Local Radiation Balance sensitivity experiment over 
the land and ocean in each hemisphere. Red: summer mean (JJA for the Northern Hemisphere; DJF for the 
Southern Hemisphere); blue: winter mean (DJF for the Northern Hemisphere; JJA for the Southern Hemi¬ 
sphere); black: annual mean. Annual mean standard deviation for each plot also shown. 


2.2.1 Overall Spatial Structure 

Figure 2.5 shows the first two global annual mean surface temperature EOF-modes for the 
Local Radiation Balance, All ON and GREB-CO 2 sensitivity experiments. In all cases, the 
dominant spatial pattern (> 88% explained variance) of the GREB model is strongly asym¬ 
metric between the hemispheres. Furthermore, the magnitude of this hemispheric pattern is 
larger in the All ON sensitivity experiment, suggesting that the GREB climate acts to en¬ 
hance this asymmetric signal. Similarly, greenhouse forcing from variable atmospheric CO 2 
is not of sufficient strength to disrupt this dominant pattern (figure 2.5e). C02-forcing does 
produce a globally synchronised signal, recorded in the monopole EOF-2 for this experiment; 
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however, the asymmetric signal remains clearly dominant in all GREB simulations. As this is 




fundamentally inconsistent with the globally-synchronised glacial signal, the GREB model 
therefore cannot reproduce the defining spatial structure of the observed Pleistocene ice ages. 
Further work is required to extract the precise mechanism through which the insolation signal 

is transformed into globally synchronised glacial cycles. 

Local Radiation Balance 
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Figure 2.5: The first two global annual mean surface temperature EOF-modes for three GREB simulations, 
top: Local Radiation Balance; middle: All ON; lower: GREB-CO 2 . 


2.2.2 Seasonal and Annual Mean Structure 

Figure 2.6 shows the spatial structure of surface temperature variability under each of the 
GREB simulations. The Local Radiation Balance climate variability is largely a function of 
latitude in annual mean and a contrasting land/sea response in the seasonal cycle. Annual 
mean variability is small and increases towards the poles. The seasonal response is much 
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more significant, with land areas varying much more than ocean areas. Atmospheric circula¬ 
tion (figure 2.6c & d) smooths this land/sea contrast in variability, particularly around coastal 
land regions. In general, this leads to an increase in temperature variability over the oceans 
and a decrease over land. 

The ice albedo feedback is most active in regions that experience seasonal ice and snow 
cover (Dommenget & Floter, 2011). In the GREB model, the strongest effects occur over the 
Southern and Arctic oceans, where there is large seasonality in sea ice cover, and the Tibetan 
Plateau in central Asia. Seasonal variability in these areas is significantly increased when the 
ice-albedo feedback is incorporated into the GREB model (figure 2.6f)- The Southern Ocean 
also experiences a small increase in annual mean temperature variability with the addition of 
the ice albedo feedback. The GREB model is most sensitive to the water vapour feedback 
in regions where atmospheric water vapour concentration is very low for at least part of 
the year (e.g. polar regions or desert areas) (Dommenget & Floter, 2011). This is because 
there is a logarithmic relationship between the water vapour content of the atmosphere and 
its emissivity. The effect of the feedback over the tropics, where the water vapour content 
of the atmosphere is high throughout the year, is therefore smaller in magnitude than over 
the high latitudes (figure 2.6g). The water vapour feedback primarily enhances the seasonal 
cycle of regions where the ice albedo feedback is active, leading to significant variability in 
the Southern Ocean and central Asia (figure 2.6h). 
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Figure 2.6: Surface temperature variability (standard deviation) for the first four GREB sensitivity experiments. 
Left column: annual mean variability; right column: seasonal cycle (June - Dec) variability. Note: colourbar 
for the seasonal cycle variability is twice as large as for annual mean. 
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Chapter 3 


OEMP & The Paleoclimate Record 


This chapter will incorporate the insights gained in chapters 1 & 2 to address the central 
question of this thesis: how does a spatially and seasonally asymmetric forcing lead to large- 
scale globally synchronised glacial cycles? In doing so, it will present and apply a simple 
insolation model based on the raw orbital parameter values. The foundation of this ‘obliquity 
and eccentricity-modulated precession’ (OEMP) model is precession forcing, the amplitude 
of which is modulated by the effects of eccentricity and obliquity. Using the OEMP model, 
it is possible to draw conclusions as to how the climate and each of the parameters interact to 
produce global annual mean surface temperature variability in the GREB model. The model 
can then be extended and compared with the observed ice age signal preserved in the 5 18 0 
record. We can therefore make interpretations as to how the Earth’s climate and Milankovitch 
forcing interact to produce globally synchronised variability over glacial-interglacial cycles. 

3.1 The OEMP Model 

The OEMP model predicts that the global mean climate evolution over the Pleistocene is 
essentially a product of orbital forcing. It is therefore assumed that a combination of each 
of the three parameters, weighted by their relative importance, will produce a time series 
comparable both to the global mean surface temperature variability of the GREB model, and 
to variations in benthic 5 18 0 in marine sediments. As each GREB sensitivity experiment, 
and the 5 ls O record, incorporates a unique combination of climate processes and feedbacks, 
the relative importance of each of the orbital parameters will similarly be unique to each ex- 
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periment. The OEMP model is therefore tuned to each climate signal, modelled or observed. 
The exact nature of this tuning can provide insight as to how insolation variability interacts 
with the climate to produce each distinct global response. The OEMP model is calculated 
with equation 3.1, which describes a simple linear combination of eccentricity (e), preces¬ 
sion (p) and obliquity (o), as defined in Huybers and Eisenman (2006). e' and o' indicate 
the anomalies of eccentricity and obliquity. Before combination, the OEMP model is tuned 
to a specific climate signal via the values taken by three scaling factors: a p , a e & a 0 . A 
graphical representation of each of the three terms in equation 3.1 for the Local Radiation 
Balance and All ON GREB sensitivity experiments is presented in figure 3.1. 


OEMP — e*a p * sin (p) + e' * a e + o' * a 0 (3.1) 


a 1 st Term: EMP 
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2 nd Term: Eccentricity Anomaly (scaled) 
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3 rd Term: Oblquity Anomaly (scaled) 



Figure 3.1: Components of the OEMP timeseries for the GREB Local Radiation Balance and All ON simula¬ 
tions. The overall insolation model (OEMP, plot d) is the sum of plots a - c. EMP (plot a) is identical for both 
GREB experiments. Note: the units for all components are dimensionless, but indicate the relative importance 
of each term. 


Eccentricity and obliquity are incorporated into the OEMP equation as anomalies as a 
consequence of the conclusions of section 1.4. The strongest interaction of both parameters 
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with precession is to increase the insolation anomaly for large values (a more elliptical orbit 
or a larger tilt) and to decrease the anomaly for low values (a more circular orbit or a smaller 
tilt). 


3.2 OEMP and GREB Climate Variability 

The OEMP model is tuned to the GREB simulations by adjusting the three scaling param¬ 
eters such that the correlation coefficient is maximised. Figure 3.2 shows that the OEMP 
model is able to reproduce almost exactly the output of the GREB model, with a correlation 
of >0.95 produced for both the Local Radiation Balance and All ON sensitivity experiments. 


OEMP vs. GREB Local Radiation Balance 



o 
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Figure 3.2: OEMP insolation model (red) compared with the GREB global annual mean surface temperature 
output (blue). Correlation (r) and values taken by the three scaling factors ( a p , a e & «„) also shown. 


The values taken by the three scaling values (shown in figure 3.2) can provide insight 
as to how the GREB model responds to each of the three orbital parameters. Firstly, a p has 
a value of +1. Physically, a sine transformation of precession (equation 3.1) means positive 
values correspond to perihelion occurring during the Southern Hemisphere summer. From 
the way the OEMP model is constructed, this means that global annual mean temperatures 
will therefore be larger than normal when perihelion occurs between December and Febru¬ 
ary. Chapter 2 discussed the importance of the Southern Hemisphere oceans in transforming 
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a seasonal precession forcing into a large annual mean signal. The value taken by a p here 
supports this conclusion, suggesting the the global mean GREB signal is led by the Southern 
Hemisphere. 

The values taken by the remaining scaling factors signal the relative importance of 
eccentricity and obliquity to the overall GREB climate response. The Local Radiation Bal¬ 
ance simulation is dominated by the eccentricity parameter, which is approximately twice 
as significant as EMP or obliquity to the overall OEMP model (see figures 3.1 & 3.2). The 
dominance of this second term in the OEMP model emerges because no climate processes 
capable of significantly enhancing the seasonally and regionally restricted forcing pattern 
of precession or obliquity are present in this GREB simulation. The absolute strength of 
eccentricity forcing is small, however, which is mirrored in the small global mean surface 
temperature anomalies that occur in this experiment (±4°C, see figure 3.2). By compari¬ 
son, the All ON simulation is less influenced by the eccentricity parameter, but the obliquity 
parameter is relatively more important. This is because internal climate feedbacks strongly 
enhance the seasonal forcing patterns of precession and obliquity. This amplification also 
means that surface temperature variability in the All ON simulation is far greater in magni¬ 
tude than in the Local Radiation Balance simulation (figure 3.2). 

The OEMP insolation model can be similarly optimised for the GREB-CO? sensitivity 
experiment (see figure 3.3). Despite incorporating an additional external forcing mechanism, 
the correlation between the OEMP model and the GREB global mean output remains high 
(>0.8). The values taken by the scaling parameters suggest that the influence of the eccen¬ 
tricity parameter is strongly enhanced in this simulation. However, as this GREB simulation 
also incorporates CO 2 forcing, this is likely a relic of the ~ 1 OOkyr cycles present in the 
atmospheric CO 2 record over the past several hundred thousand years (Liithi et al., 2008). 

The correlation between the OEMP model and the GREB climate increases signifi- 
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OEMP vs. GREB-C0 2 



Figure 3.3: OEMP insolation model (red) compared with the GREB CCF-Forced global annual mean surface 
temperature output (blue). Correlation (r) and values taken by the three scaling factors (a p , a e & a 0 ) also 
shown. 

cantly if the ~413kyr eccentricity signal is removed from the OEMP calculation. Variability 
on this frequency has been pointed out as notably absent from the observed paleoclimate 
record (Berger, 1978; Ruddiman, 2003; Willeit et al., 2019). Several physical explanations 
for this missing signal have been proposed (e.g. Ruddiman, 2003), however this charac¬ 
teristic of glacial cycles remains poorly understood. Regardless, the results presented here 
support the theory that the ~413kyr eccentricity frequency is of relatively little significance 
in contributing to the global annual mean climate evolution across glacial-interglacial cycles. 

3.3 OEMP and Observed Pleistocene Climate Variability 

The OEMP model is highly effective in reproducing the global annual mean surface temper¬ 
ature output of the GREB model. However, several features of the Earth’s climate system 
relevant on glacial timescales are not simulated in the GREB model. Firstly, the GREB 
ocean extends to a maximum depth of three times the seasonally-shifting mixed layer depth, 
with deeper ocean layers not simulated. Similarly, the GREB model does not include an ice 
sheet model, nor a fully interactive carbon cycle. These additional processes and feedbacks 
are thought to play a significant role in the amplitude and time-scale behaviour of glacial- 
interglacial cycles (Abe-Ouchi et al., 2013; Ellis & Palmer, 2016; Berger, Li, & Loutre, 
1999). Despite these limitations, the OEMP model is again able to successfully reproduce 
the major features of the 5 18 0 signal, with a correlation coefficient of -0.723. As for the 
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GREB-CO 2 simulation, removing the ~413kyr eccentricity signal significantly improves 


the overall correlation between the two timeseries (see figure 3.4). 


OEMP vs. 4 18 0 



OEMP (with 413kyr Ecc. Cycle) vs. 4 18 0 



Figure 3.4: OEMP insolation model (red) compared with the f) lx O stack (blue). Upper: 413kyr eccentricity 
cycle removed from the OEMP model; lower: 413kyr eccentricity cycle preserved. Correlation (r) and values 
taken by the three scaling factors ( a p , a e & <X 0 ) also shown. 

The scaling parameters demonstrate that the Earth’s climate responds to orbital forcing 
in a fundamentally different way to the GREB model. Most significantly, a p takes a value 
of -1. Under this transformation, the global mean temperature anomaly is now larger when 
perihelion occurs during the Northern Hemisphere summer, meaning that the global mean 
climate response to insolation forcing is led by this Hemisphere. This is the opposite scenario 
to that of the GREB model, but agrees with the predictions of the prevailing scholarly opin¬ 
ion (Elkibbi & Rial, 2001). Current theory holds that the margins of continental Northern 
Hemisphere ice sheets are positioned at approximately 65°N; thus, these ice sheets are most 
sensitive to summer insolation variability at this latitude. The results of the OEMP model 
support this theory, suggesting that insolation forcing over the Northern Hemisphere drives 
the global response. As the GREB model does not simulate ice sheet dynamics, the opposite 
scenario develops and it is the Southern Hemisphere that ultimately determines the GREB 
climate. By comparison, the direct effect of eccentricity is identical to that of the prescribed 
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CO 2 experiment, while the effect of obliquity is significantly enhanced. This is because the 
relative importance of the obliquity parameter is heavily reliant on internal feedback mech¬ 
anisms. In addition to the ice albedo and water vapour feedbacks of the GREB model, the 
5 18 0 record is strongly influenced by feedbacks in the carbon cycle and through ice sheet 
dynamics. The spatial distribution of ice sheets in the very high latitudes in particular render 
them highly responsive to obliquity forcing. The net effect is that this parameter is almost 
twice as significant to the time evolution of the 5 18 0 timeseries than to the GREB model. 

The success of the OEMP model in reproducing the major features of the <5 I8 0 record 
suggests that orbital variability is the primary forcing mechanism and ultimate pacemaker 
of the glacial-interglacial cycles characteristic of last 2.5myr. The model itself consists only 
of a scaling and linear combination of the raw parameter data with a very small set of as¬ 
sumptions as to how the climate should interact with this forcing. This implies that internal 
climate processes principally act to (1) transform a hemispheric insolation forcing into an 
global mean signal; and to (2) amplify this underlying forcing through positive feedbacks. 
In relation to the observed ice age signal in the 5 18 record, ice sheet dynamics and the carbon 
cycle likely perform important functions to this effect. 
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Conclusion 


The results of this thesis demonstrate that the dominant structure of orbital insolation vari¬ 
ability is asymmetric between the hemispheres. This signal is primarily produced by varia¬ 
tions in precession, modulated by the eccentricity of the Earth’s orbit. The direct influence 
of the eccentricity and obliquity parameters act to spatially and seasonally rearrange this 
dominant precession signal, particularly in the high latitudes. The response of the GREB 
model to orbital insolation forcing mirrors this asymmetrical structure, with a principal sur¬ 
face temperature response of opposite sign between the hemispheres. This spatial structure 
persists even with the inclusion of prescribed variability in atmospheric CO 2 concentrations. 
The GREB model is therefore unable to successfully reproduce the globally synchronised 
glacial-interglacial cycles recorded in paleoclimate data sources. However, the deconstruc¬ 
tion of the GREB climate response indicates that the oceans, atmospheric circulation, the 
ice-albedo feedback and the water vapour feedback each significantly influence the regional, 
seasonal and annual mean structure of temperature variability on glacial timescales. Fur¬ 
thermore, the OEMP model demonstrates the dominant role of the Northern Hemisphere in 
driving variations in benthic 5 18 0 isotope ratios over the past lmyr. This is opposite to the 
surface temperature signal in the GREB simulations, suggesting that climate feedbacks that 
are not simulated in the model play an important role in determining the hemispheric re¬ 
sponse to insolation forcing. Our results suggest that these feedbacks must be characterised 
by a large inertia that allows for a strongly seasonal forcing to be transformed into a coherent 
annual mean signal. The dynamic response of ice sheets and the global carbon cycle con¬ 
stitute plausible physical mechanisms that satisfy this condition. Future work should focus 
on determining the precise Earth system mechanisms through which a globally synchronised 
climate signal can be produced in response to a fundamentally asymmetric forcing pattern. 
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Appendix 1: Literature Review 


Introduction 

First described mathematically by Milankovitch (1941), the astronomical theory of climate 
states that large-scale glacial fluctuations are controlled by secular variations in the Earth’s 
orbit around the sun. This variability is generated by long-term oscillations in three key 
orbital parameters, each of which affects the regional and spatial distribution of insolation 
across the Earth’s surface. The first of these, eccentricity, measures the elliptical shape of 
the Earth’s orbital path and varies at periods of 413kyr (413,000 years), 125kyr and 95kyr. 
The obliquity parameter, by comparison, oscillates with a 41kyr cycle and describes the an¬ 
gle between the Earth’s axis of rotation and the plane of its orbit. Finally, the precession of 
the equinoxes dictates how the seasons are positioned in relation to perihelion, fluctuating 
with key periodicities of 23kyr and 19kyr. The spectra of the orbital parameters reveal cycli¬ 
cal behaviour on multiple frequencies, resulting in a complex insolation signature spanning 
geological time scales. Hays, Imbrie, and Shackleton (1976) published the first concrete 
evidence that the Earth’s climate is linked to orbital variability, showing that the spectra of 
<5 18 0 isotopic ratios from marine sediment cores yield peaks at frequencies that closely align 
with those of the three parameters. Marine 5 18 0 ratios are interpreted as proxies for global 
ice sheet extent, as ice sheets are enriched in the lighter 160 oxygen isotope, which forms a 
larger fraction of the snow that falls over the higher latitudes. The prevailing interpretation 
of this spectral alignment is that orbital variability is the ultimate source of continental-sc ale 
glacial cycles (Shackleton, 2000). 

Despite the broad scientific consensus as to the validity of Milankovitch theory, the 
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physical connection between orbital forcing and the Earth’s climate remains poorly under¬ 
stood. Modelling initiatives have been successful in simulating several key features of the 
ice ages (e.g. Pollard, 1983; Paillard, 1998; Tziperman, Raymo, Huybers, & Wunsch, 2006; 
Huybers, 2011; Parrenin & Paillard, 2012). No models to date, however, have been able 
to simulate all characteristics of the observed paleoclimate data, nor has a conclusive set of 
physical processes or feedbacks been established that is able to fully explain all aspects of the 
ice age climate signal (Willeit, Ganopolski, Calov, & Brovkin, 2019). This review will high¬ 
light the key areas of uncertainty that remain. The challenges are divided into three key areas: 
(1) the structure of orbital insolation variability; (2) the non-linear climate response; and (3) 
the transition from a dominant period of 41kyr to lOOkyr roughly one million years ago (the 
Mid-Pleistocene Transition or MPT). This review concludes that although the Earth’s cli¬ 
mate is clearly linked to long-term orbital variability, the precise relationship between the 
pattern of insolation forcing and the observed climate response requires further explanation. 

Orbital Insolation Variability 

Each of the three orbital parameters oscillates on multiple frequencies and influences in¬ 
solation on different regional and spatial scales (see figure 1). As a result, there remains 
significant uncertainty as to which parameter(s) is the primary driver behind the ice age 
signal preserved in the paleoclimate record. Numerous models have been shown to effec¬ 
tively simulate the key features of glacial cycles, but these lack agreement as to whether 
eccentricity (e.g. Rial, 1999; Ruddiman, 2003), obliquity (e.g. Huybers & Wunsch, 2005), 
precession (e.g. Ridgwell, Watson, & Raymo, 1999; Ellis & Palmer, 2016), or some combi¬ 
nation thereof (e.g. Ganopolski & Calov, 2011; Huybers, 2011; Parrenin & Paillard, 2012) 
can be considered the dominant forcing. For this reason, our research will investigate how 
each of the three orbital parameters, both individually and collectively, influence the regional 
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and seasonal distribution of insolation on millennial timescales. 



Figure 1: Insolation variability over the past 1 million years, a mean July insolation at 65°N at the top of the 
atmosphere (W/m 2 ); b eccentricity (dimensionless); c obliquity (degrees); and d climatic precession (dimen¬ 
sionless). From (Berger et al., 1999) 


Of the three parameters, only eccentricity is able to influence the total amount of ra¬ 
diation that reaches the top of the atmosphere each year (Berger, 1978). As well as this, the 
~lOOkyr ice age cycles of the most recent 800kyr are seemingly aligned with one of the two 
main eccentricity periodicities, suggesting that this parameter has a significant influence over 
the climate. However, the direct effect of eccentricity on insolation at lOOkyr timescales is 
far too small (~ 1 W/m 2 or < 1 % of total insolation variability) for it alone to be the key driver 
of glaciation events (Berger, 1978; Elkibbi & Rial, 2001). Similarly, eccentricity varies at 
frequencies of 413kyr, 125kyr and 95kyr, none of which coincides exactly with the power 
spectra of marine 5 18 0 (see figure 2) (Imbrie et al., 1993). Various studies have shown that 
any particular spectral peak does not necessarily signify the presence of such a periodicity in 
the climate (Maslin & Brierley, 2015; Berger, Melice, & Loutre, 2005). The lOOkyr peak in 
the <5 18 0 spectrum could be the result of overlapping cycles of precession (Ridgwell et al., 
1999) or obliquity (Huybers & Wunsch, 2005), and thus can be better interpreted as a relic 
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of the statistical technique itself (Maslin & Brierley, 2015). 



Figure 2: Benthic 5 18 0 from an Atlantic (Imb84) and Pacific (OJsox96) marine sediment core, a time evolution 
of 5 18 0 over the most recent 700kyr. b the power spectra of a, with 400kyr, 125kyr, 95kyr (eccentricity); 41kyr 
(obliquity) and 23kyr and 19kyr (precession) periodicities indicated. From (Maslin & Brierley, 2015) 

Despite having no influence over the total annual insolation budget of the planet, pre¬ 
cession and obliquity dominate total insolation variability due to orbital forcing. These pa¬ 
rameters achieve this by controlling how insolation is distributed across latitudes and sea¬ 
sons. Precession acts to enhance the seasonal cycle in one hemisphere (colder winters and 
warmer summers) while simultaneously dampening the seasonal cycle in the other (warmer 
winters and cooler summers). As either side of the planet experiences a forcing anomaly 
of opposite sign, precessional forcing is said to be ‘out of phase’ across the hemispheres 
(Raymo & Huybers, 2008). Unlike precession, obliquity is in-phase between the two hemi¬ 
spheres, meaning an increase in obliquity increases the seasonal gradient in both hemispheres 
simultaneously. The largest effects of obliquity are felt in the high latitudes where the sea¬ 
sonal gradient is most pronounced (Berger 1978). Thus, precession dominates low latitude 
insolation variability whereas obliquity dominates near the poles, with the two parameters 
roughly equal in strength at around 65°north and south of the equator (Kukla & Kukla, 1972). 

Insolation Metrics 

In the context of Pleistocene glacial cycles, it is important to determine the exact latitude 
and season at which the climate is most sensitive to insolation variability. Milankovitch 
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(1941), building on the work of Koppen and Wegener (1924), argued that summer insolation 
in the high northern latitudes (usually defined as June 21 at 65°N) ultimately forces ice sheet 
variability. This theory states that the mass balance of an ice sheet is largely a function of 
the summer ablation rate at its margins, which are estimated to be located at approximately 
65°N (Oerlemans, 1991). This definition has become the most widely accepted insolation 
metric used to analyse the relationship between orbital forcing and glacial cycles (Weertman, 
1976; Oerlemans, 1991). However, the effectiveness of this metric has not been adequately 
established (see figure 3), and several alternatives have been proposed. Kukla (1975) argues 
that ice sheets are more sensitive to insolation changes in autumn than summer. Alternatively, 
(Huybers, 2006) argues that the integration of insolation over an entire summer provides 
a better measure of orbitally-induced ablation rates. A variation of this metric is ‘caloric 
summer half-year insolation’, defined as the period during which any day receives more 
insolation than any day in the winter half-year. This accounts for the fact that summers that 
fall near perihelion (due to changes in precession) are also shorter than usual due to the Earth 
travelling faster through that section of its orbit (Tzedakis, Crucifix, Mitsui, & Wolff, 2017). 
Which of these insolation metrics (or some other alternative) is most effective for ice age 
climate analysis remains a key area of uncertainty. 

The Climate Response 

Figure 3 compares the variations in June 21 daily mean insolation at 65°N with a 5 18 0 iso¬ 
topic ratio extracted from a Pacific sediment core. This relationship is regularly presented 
as clear evidence of a connection between orbital forcing and global ice volume (Berger et 
al., 1999; Ganopolski & Calov, 2011; Ellis & Palmer, 2016). It is argued that interglacial 
periods (labelled with arrows and defined as the periods at which the benthic 5 18 0 crosses 
a 3.68%c threshold) occur concurrently with summer 65°N insolation maxima. Immediately 
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upon inspection, however, several inconsistencies in this relationship become apparent. Al¬ 
though a handful of interglacials coincide roughly with large insolation maxima (blue shaded 
regions), there are just as many interglacials that are associated with only moderate insola¬ 
tion peaks (red shaded regions). Similarly, the lOOkyr cycles argued to dominate this period 
lack any consistency, with clear and sustained deviations from this apparent regime (e.g. be¬ 
tween 700kyr and 500kyr BP; orange arrow). Based on total 21 June insolation variability 
alone, it is impossible to predict the time evolution of benthic 5 18 0 (Tzedakis et al., 2012). 
As this metric forms the basis of the great majority of research and modelling into Pleis¬ 
tocene glacial cycles, a more robust analysis of the regional and seasonal structure of orbital 
insolation forcing is required. This will be a major focus of my honours research. 


No lOOkyr period 

◄-► 



Age (kyr bp) 

Figure 3: a benthic 8 18 0 stack for the past 800kyr. b Daily mean insolation at 65°N on June 21. Marine Isotope 
sub-Stages displayed above a, arrows show interglacials, defined as the point at which the benthic <5 18 0 crosses 
a 3.68%o threshold. From (Tzedakis et al., 2017) 

Irrespective of the insolation metric chosen by researchers, the irregular time evolution 
of global ice volume over the past 800kyr demonstrates that the Earth’s climate-cryosphere 
interacts with the insolation forcing in a strongly non-linear and, as yet, poorly understood 
manner. The only significant orbital insolation is regional and seasonal, a signal inconsistent 
with the coherent global response recorded in the paleoclimate record (Doughty et al., 2015). 
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During the last ice age, for example, both northern and southern hemisphere ice sheets have 
been shown to have varied in concert (Mercer, 1984; Denton et al., 1999). As precessional 
insolation forcing is out of phase between the hemispheres, this implies that one hemisphere 
does not respond to regional insolation forcing in the manner predicted by Milankovitch the¬ 
ory. These inconsistencies between the orbital forcing and the climate response indicate that 
some physical process(es) internal to the climate system must act to transform a spatially and 
temporally restricted insolation pattern into global glacial cycles. Several feedback mech¬ 
anisms have been proposed as potential solutions to this problem. These have been tested 
in various models, and some are effective in reproducing the major characteristics of the ice 
age signal, including the ‘sawtooth’ pattern of time evolution and approximate lOOkyr pe¬ 
riodicity (Oerlemans, 1980; Pollard, 1982; Berger et al., 1999). In spite of these advances, 
however, the precise processes responsible for the non-linear climate response to orbital 
forcing remain elusive. 

Atmospheric CCb concentrations, effectively uniform across the globe, provide a clear 
mechanism for transforming a regional insolation signal into a global glaciation event. How¬ 
ever, the role of greenhouse gases, and carbon dioxide specifically, in contributing to Pleis¬ 
tocene climate variability is poorly understood. Almost all studies acknowledge that CO 2 
likely plays a role in amplifying the ice age signal; however, models are split between as¬ 
serting that it is either critical (Crowley & Hyde, 2008), somewhat critical (e.g. Ganopolski 
& Calov, 2011) or non-determinative (e.g. Abe-Ouchi et al., 2013) in the time evolution of 
the glacial cycles. While atmospheric CO 2 can provide insight as to the amplification of ice 
ages into large-scale global events, it cannot explain the asymmetrical time evolution of the 
glacial cycles. An enduring physical hypothesis for the non-linear climate response points to 
the delayed response of the lithosphere to changes in ice sheet mass (Ganopolski & Brovkin, 
2017). This regolith-ice sheet interaction works to keep ice sheets at higher, colder altitudes 
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as they accumulate mass, and at lower altitudes as they ablate. Models have shown that this 
feedback can reproduce the observed ‘sawtooth’ pattern of sustained accumulation prior to 
a glacial maximum, and rapid ablation at a glacial termination (Abe-Ouchi et al., 2013). An 
alternative physical explanation relates to a dust-albedo feedback. This theory states that 
glacial maxima are characterised by increased aridity and soil erosion rates, which com¬ 
bine to produce large dust-storms that deposit particulate matter onto continental ice sheets 
(Ganopolski & Calov, 2011; Ellis & Palmer, 2016). This lowers both the albedo of an ice 
sheet and the insolation required to trigger its rapid ablation and initiate an interglacial pe¬ 
riod. Numerous additional feedback mechanisms have been put forward as essential to the 
non-linear response of the climate to astronomical forcing, most notably ice-sheet calving 
(Pollard, 1982) and ocean dynamics (Deblonde, Peltier, & Hyde, 1992). My research will 
attempt to shed light on this large area of uncertainty by simulating the global climate over 
the entire Pleistocene using a simple globally-resolved energy balance (GREB) model. With 
the GREB model, it is possible to simulate the climate response to external forcing while 
removing the effects of key climate feedbacks. Similarly, the GREB model can be run with 
either variable or constant atmospheric CO 2 concentrations to establish its role as either a 
driver or amplifier of the ice age signal. These sensitivity experiments will provide insight 
as to how the orbital forcing signal is scaled and transformed by internal climate processes. 

The Mid-Pleistocene Transition 

Cycles of large-scale glacial expansion and retreat have occurred since the beginning of the 
Pleistocene, roughly 2.6 million years ago. Prior to this, it is theorised that elevated atmo¬ 
spheric CO 2 concentrations suppressed large-scale Northern Hemisphere ice sheets (Lunt, 
Foster, Haywood, & Stone, 2008). The first ice age cycles, recorded by ice-rafted debris in 
the North Atlantic, were smaller in magnitude that those of the late Pleistocene and charac- 
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terised by a regular 41kyr periodicity (Raymo & Nisancioglu, 2003). Between ~1.25myr 
and ~0.8myr ago, however, these cycles abruptly shifted to an irregular lOOkyr periodicity 
during the Mid-Pleistocene transition (see figure 4). 



Figure 4: A benthic 5 18 0 stack for the past 3myr, showing the transition from a dominant periodicity of 
41kyr to lOOkyr during the MPT, ~ 1 million years ago. Modified from Raymo and Huybers (2008), data from 
Lisiecki and Raymo (2005). 

Various hypotheses have been presented to account for the MPT, including changes to 
atmospheric greenhouse gas concentrations, regolith and climate-dust feedbacks, and inter¬ 
pretations of the proxy data itself. Raymo, Lisiecki, and Nisancioglu (2006) argue that the 
apparent 41kyr cycles of the early Pleistocene are largely a relic of the globally integrated 
marine 5 18 0 proxy used to estimate ice volume over the last 3-5myr. Models have shown 
that under certain conditions, enhanced ablation from precessional forcing in one hemisphere 
could be balanced by equal accumulation in the other (Raymo et al., 2006). This would ef¬ 
fectively remove the precessional signal from a globally integrated proxy record, leaving 
only the 41kyr periodicity of obliquity. However, this hypothesis does not explain the abrupt 
end to this early glacial regime and the inception of the ‘lOOkyr world’. 

Numerical solutions to the variations in the three Milankovitch parameters have shown 
that orbital insolation forcing has not changed significantly for millions of years (Berger, 
1978; Huybers & Eisenman, 2006), implying that the MPT must have been forced by a 
shift in internal climate processes (Ellis & Palmer, 2016). Clark and Pollard (1998) point to 
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the gradual process of regolith removal over the Northern Hemispheric continents as a key 
driver of the MPT. The authors argue that thick sediments facilitate increased rates of basal 
sliding, leading to thinner ice sheets with larger ice losses during periods of summer ablation. 
(Ganopolski & Calov, 2011) test this hypothesis with a coupled climate-ice sheet model 
of intermediate complexity, finding that the lOOkyr glacial cycle is removed when thick 
sediments are prescribed across all land areas. An alternative hypothesis argues that a similar 
threshold effect, driven primarily by steady reductions in atmospheric CO 2 concentrations, 
ultimately brought the dominance of the 41kyr cycle to an end (Berger et al., 1999; Elkibbi 
& Rial, 2001; Ruddiman, 2003). This theory argues that a gradual decrease in atmospheric 
CO 2 driven by tectonic and oceanic processes was accompanied by global cooling, which 
eventually allowed large ice sheets to develop and persist even during periods of relatively 
high insolation (Willeit, Ganopolski, Calov, Robinson, & Maslin, 2015). Although high- 
resolution CO 2 data does not extend this far into the Earth’s past, the trend observable in 
marine 5 18 0 ratios is suggestive of a gradual cooling throughout the Pleistocene (Ruddiman, 
2003). 

Conclusion and Research Directions 

Despite decades of research and numerous modelling initiatives, the exact relationship be¬ 
tween orbital forcing and the glacial-interglacial cycles of the past ~2.6 million years re¬ 
mains largely uncertain. Specifically, the regional, seasonal and long-term structure of or¬ 
bital insolation variability over glacial timescales is poorly defined. Similarly, there is little 
agreement as to how the Earth’s climate interacts with this insolation forcing to produce the 
observed ice age signal, in particular its irregular time evolution and the regime shift of the 
Mid-Pleistocene Transition. My research will focus on these key areas. Firstly, by offering 
a statistical analysis of the seasonal, spatial and time-scale behaviour of the Milankovitch 
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forcing; and secondly, by analysing the climate response of the GREB model to this orbital 
forcing. Due to the simplicity of the GREB model, it is possible to run multiple simulations 
over the entire Pleistocene. Simulations of this length are currently lacking in the literature, 
and will provide insight as to how the climate responds to orbital forcing on the scale of 
millions of years. This is particularly relevant in confirming whether the MPT is the result 
of internal climate processes, or simply a response to long-term fluctuations in orbital vari¬ 
ability. These experiments will provide a more complete understanding of the spatial and 
temporal structure of the orbital forcing, as well as how it relates to climate variability on 
multiple timescales. 
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